function X = solve_beta(beta,gamma,psi,mu_C,sig_C,P,Pnodis)


alpha   = 1-gamma;
rho     = 1-1/psi;
theta   = alpha/rho;

nS      = length(P);
stat    = Pnodis^1e6;
stat    = stat(1,:)';

% Solve value function
lnWC        = ones(nS,1);
opt         = optimset('Display','off','TolFun',1e-12,'MaxFunEval',1e6);
f_lnWC      = @(lnWC) theta*log(beta) + alpha*mu_C + 1/2*alpha^2*sig_C.^2 + log(P*(exp(lnWC)+1).^theta) - theta*lnWC;
lnWC        = fsolve(f_lnWC,lnWC,opt);
lam         = exp(lnWC);


% Pre-compute some SDF related objects
S           = ((lam'+1)./lam).^(theta-1);


% Risk-free rate
A           = beta^theta * bsxfun(@times,S,exp(-gamma*mu_C+0.5*gamma^2*sig_C.^2));
bond_f      = NaN(nS,120);
bond_1      = ones(nS,1);
for i=1:120
    bond_1      = (P.*A)*bond_1;
    bond_f(:,i) = bond_1;
end
Rf      = 1./bond_f; % [nS,120]

%Rf12 = (stat'*Rf([1,3],1) - 1)*12;
Rf12 = stat'*Rf([1,2],12) - 1;

X = Rf12 - .01;